Morphological stability of electromigration-driven vacancy islands 
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The electromigration-induced shape evolution of two-dimensional vacancy islands on a crystal 
surface is studied using a continuum approach. We consider the regime where mass transport is 
restricted to terrace diffusion in the interior of the island. In the limit of fast attachment/detachment 
kinetics a circle translating at constant velocity is a stationary solution of the problem. In contrast 
to earlier work [O. Pierre-Louis and T.L. Einstein, Phys. Rev. B 62, 13697 (2000)] we show 
that the circular solution remains linearly stable for arbitrarily large driving forces. The numerical 
solution of the full nonlinear problem nevertheless reveals a fingering instability at the trailing end 
of the island, which develops from finite amplitude perturbations and eventually leads to pinch-off. 
Relaxing the condition of instantaneous attachment/detachment kinetics, we obtain non-circular 
elongated stationary shapes in an analytic approximation which compares favorably to the full 
numerical solution. 

PACS numbers: 05.45.-a, 68.65.-k, 66.30.Qa, 68. 35. Fx 
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I. INTRODUCTION 



Much of the diversity of natural shapes in the inan- 
imate world is the result of morphological instabilities. 
The paradigmatic example is the Mullins-Sekerka insta- 
bility, in which a spherical solid nucleus in an undcrcooled 
melt forms lobes and petals which eventually develop into 
a delicate dendritic pattern and its two-dimensional 
analogs that can be observed in the growth of islands 
on crystal surfaces [2|. These systems share the com- 
mon mathematical structure of moving boundary value 
problems, in which an interface evolves in response to 
the gradient of some continuous field defined in the spa- 
tial domains that it separates. In the context of two- 
dimensional crystal surfaces, the interfaces are atomic 
height steps separating different terraces, and their mo- 
tion is governed by the attachment and detachment of 
the adsorbed atoms (adatoms) [!, 0, H[ . 

A rich variety of two-dimensional morphological in- 
stabilities has been observed on the surfaces of current- 
carrying crystals, where an electromigration force induces 
a directed motion of adatoms 0, 0] ■ The microscopic ori- 
gin of this force is a combination of momentum transfer 
from the conduction electrons (the 'wind force') and a 
direct effect of the local electric field [1] . On stepped sur- 
faces vicinal to Si(lll), electromigration has been found 
to cause step bunching 0, step meandering [Tcj | . step 
bending 11 1 and step pairing [12| instabilities. In addi- 
tion, single layer adatom islands on silicon surfaces have 
been seen to drift under the influence of electromigration 
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In the present paper we focus on the morphological 
stability of single layer vacancy islands driven by elec- 
tromigration. We build on the work of Pierre-Louis 
and Einstein (PLE) [III, who introduced a class of con- 
tinuum models for island electromigration. The dif- 
ferent models are distinguished according to the domi- 
nant mechanism of mass transport, which can be due to 
periphery diffusion (PD) along the edge of the island, 
terrace diffusion (TD), or two-dimensional evaporation- 
condensation (EC), i.e. attachment-detachment, kinet- 
ics. In the PD regime the dynamics of the island edge 
is local, while in the TD and EC regimes it is coupled 
to the adatom concentration on the terrace. In the TD 
(EC) regime diffusion is slow (fast) compared to the at- 
tachment/detachment processes, as reflected in the mag- 
nitude of the kinetic lengths 



d± = D/k± 



(1) 
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defined as the ratio of the surface diffusion constant D to 
the rates of adatom attachment to a step from the lower 
(fe+) or upper (fc_) terrace, respectively. The two rates 
generally differ, because step edge barriers suppressing 
attachment across descending steps are ubiquitous on 
many surfaces, leading to k+ > k- 0- 

As a common feature of the models of PLE, the elec- 
tromigration force is taken to be of constant direction 
and magnitude everywhere, which implies in particular 
that it is not affected by the presence and the shape of 
the island. This is motivated by the fact that the island 
constitutes a small perturbation in the morphology of the 
crystal, which is not expected to substantially change the 
distribution of the electrical current in the bulk. Detailed 
atomistic calculations do in fact show that the electromi- 
gration force is modified in the vicinity of a step 0, [ljj , 
which may also be incorporated into a continuum model 
[l8l |. but this is a higher order effect that can be ne- 
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glected on the present level of description. The situation 
is completely different for electromigration-driven macro- 
scopic voids in metallic thin films, which can be modeled 
using a closely related two-dimensional continuum ap- 
proach 0, [H, [M H H ll, [H . Since the void inter- 
rupts the current flow, the effect of the void shape on the 
current distribution is an essential part of the analysis. 
In the following we refer to this problem as void migra- 
tion, to be distinguished from island migration under a 
constant force. 

Apart from the work of PLE, analytical results con- 
cerning the morphological stability of electromigration- 
driven two-dimensional shapes have been obtained only 
in the PD regime. In the absence of crystal anisotropy the 
basic solution is then a circle moving at constant veloc- 
ity [19(. In the case of island electromigration the circle 
becomes linearly unstable at a critical radius or critical 
driving force [20j . Beyond the linear instability station- 
ary sha pes that are elongated in the current direction 
appear |26l . |27| . Remarkably, the stability scenario for 
void migration is completely different. Voids are linearly 
stable at any size [2(J, |2lj , but they become nonlinearly 
unstable beyond a finite threshold perturbation strength, 
which decreases with increasing radius or driving force 
|22l l24l ]. Unstable voids break up into smaller circular 
voids, and non-circular stationary shapes do not exist 
(25)]. The increasing sensitivity to finite amplitude per- 
turbations can be linked to the increasing non-normality 
of the linear eigenvalue problem, which leads to transient 
growth of linear perturbations [22j • This route to nonlin- 
ear instability has been previously described for linearly 
stable hydrodynamic flows [H, [2^, [3(3], and it will be 
further discussed below in Section ITVl 
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FIG. 1: Schematic of the interior model. Adatoms detach 
from the inner boundary of the vacancy island and diffuse 
subject to the electromigration force F directed to the right. 
As a consequence, the entire island drifts to the left at speed 
V. 

In the present paper we focus on the "interior model" 
introduced by PLE. Referring to Fig. [TJ we consider a 
single vacancy island which is isolated from the surround- 
ing upper terrace by a strong step edge barrier that pre- 
vents adatoms from entering across the descending step 
(fc_ =0, d_ = 00). Diffusive motion of vacancy islands 
mediated by internal terrace diffusion has been observed 
experimentally on the Ag(110) surface [3l|. The math- 



ematically equivalent process of internal diffusion of va- 
cancies also pla ys an important role in the motion of 
adatom islands [32]. We will use the terminology ap- 
propriate for adatom diffusion inside a vacancy island 
throughout the paper. 

The mathematical description of the interior model 
leads to a moving boundary value problem on a finite 
domain, which we formulate in the next section. A key 
ingredient of our analytic work is a separation ansatz for 
the adatom concentration, which allows us to determine 
stationary island shapes and investigate their linear sta- 
bility in a simpler and more transparent way than in pre- 
vious work [la ] . The analytic approach is complemented 
by numerical simulations of the full nonlinear and non- 
local dynamics. Specifically, in Section IlIII wc compute 
non-circular stationary island shapes perturbatively in 
the parameter <5 = <i + /i?o, where Rq is the radius of 
the circle which solves the stationary problem in the TD 
limit (d+ = 0). Section [TVl is devoted to the linear sta- 
bility analysis of the circular solution for d+ = 0. In 
agreement with PLE, we find that the eigenvalues of the 
linearized problem depend only on the ratio z = i?o/£, 
where 



is the characteristic length scale associated with the elec- 
tromigration force F. However, in contrast to PLE, who 
argued (on the basis of a less extensive analysis) that the 
circle becomes unstable for z > 0.1, we show that it in 
fact remains linearly stable for all values of z. Simula- 
tions of the full nonlinear evolution in SectionlVlneverthe- 
less reveal an instability under finite amplitude perturba- 
tions, in which a finger develops at the trailing end of the 
island and eventually leads to a pinch-off. In Section IVII 
we summarize our results and discuss their significance in 
the broader context of morphological stability in moving 
boundary value problems. 

II. MODEL 

Since we assume that the mass transport on the surface 
is dominated by terrace diffusion, the main dynamical 
quantity of interest is the adatom concentration c(x, y, t) 
on the terrace. By mass conservation its time evolution 
is governed by 

<9 t c + V-j = (3) 
D 

j = -DVc + — xc, (4) 

where the current j takes into account the contributions 
from diffusion and electromigration. Since we assume the 
force to be constant, electromigration appears as a drift 
in the direction of the electric current, denoted by the 
unit vector x. 

The coupling to the periphery of the island is given 
by the boundary conditions at the step edge. Let the 
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subscripts +, — denote quantities at the lower and up- 
per terrace, respectively, n the normal pointing from the 
upper to the lower terrace, and v the normal velocity of 
the island boundary. The fluxes 



3± := TO'd 



c±v) 



(5) 



from the lower (+) and the upper (— ) terrace, respec- 
tively, towards the step are then assumed to be propor- 
tional to the deviation from equilibrium i.e., 



i± = k± (c± 



(6) 



with fc + , k- denoting the attachment rates from the 
lower and upper terrace, respectively. Here the equi- 
librium density c cq is given by the linearized Gibbs- 
Thomson relation 



(1+Tk) 



r = a^/keT, 



(7) 



where 7 denotes the (isotropic) step stiffness, a the lattice 
constant, and k the curvature of the terrace boundary. 
The validity of ([7]) requires the capillary length V to be 
small compared to the radius of curvature of the island 
boundary. Defining the kinetic lengths by |T]), we note 
that in the terrace diffusion limit (TD), where the at- 
tachment/detachment becomes instantaneous (k± — > 00, 
d± — > 0), the boundary conditions ([6]) reduce to 



(8) 



Finally, by mass conservation, the normal velocity v of 
the boundary is 



(9) 



We neglect periphery diffusion, since we are concerned 
here with the kinetic regime where terrace diffusion is the 
dominant mass transport mechanism. For the interior 
model fc_ = 0, which implies j'_ = 0, and ([H]) applies on 
the interior (lower) terrace in the TD limit. 

In the following, it is assumed that the electric cur- 
rent is in the cc-direction, i.e. x — (1,0). Moreover, 
for analytic calculations the quasistatic approximation is 
frequently used, which amounts to setting dtc = in the 
diffusion equation ©, which then reduces to 



Ac-C 1 d x c = 



(10) 



and omitting the term proportional to v in ([5]), which 
arises from the sweeping of adatoms by the moving step. 
For our purposes, the general solution of (fTU|) is most 
conveniently expressed in polar coordinates. To arrive 
at a suitable representation, we first eliminate the drift 
term breaking the rotational symmetry of (|10p via the 
ansatz 

c(x,y) = exp(— )/(— , — ), 

which leads to the Helmholtz equation Af = f. Sepa- 
ration of the latter equation yields a harmonic angular 



dependence and a modified Bessel function of imaginary 
argument I n for the radial part. The general solution of 
(fTUl) is then a superposition of the form 



c(r,9) = exp(^cos0)/(r,0), 

OO 

f(r,6)= ^2 c 7l 7 n (^)exp(m6 



(11) 
(12) 



n— — 00 



where the unknown coefficients {c n } are to be determined 
by the boundary conditions 

We further note that, in the quasistatic approximation, 
the total area A of the island is strictly conserved by the 
dynamics. For the interior model one computes 



-A 



v ds 



mi 



.)+ 



ft ds 



00 



-a- 1 V • j+ dA = 0, 
In 



using the divergence theorem, where denotes the in- 
terior domain and dfl its boundary. The last integral 
vanishes in the quasistatic approximation. This means 
that the mass exchange between the island boundary (the 
bulk) and the adatom concentration inside is always bal- 
anced. In this sense the diffusion field merely mediates 
the mass transport from one part of the boundary to 
another. 



III. STEADY STATES 

For the interior model in the TD limit (d+ = 0, cL = 
00), circular islands are steady states. More precisely, an 
island with radius Rq, constant adatom concentration 



c° (1 



c - -eq 

and drifting with constant velocity 



r ' 
Ro ■ 



V 



D 



a 2 c 



£ 1 



a 2 Co' 



(13) 



is a solution of ([I])-©; the factor (1 — a 2 c ) _1 is a cor- 
rection to the quasistatic approximation, which requires 
a 2 co <C 1. However, if the attachment is not instanta- 
neous (d+ > 0) , the circle is no longer stationary. As was 
shown in [15| |. an expansion of the interior model to sec- 
ond order in z = i?o/£ leads to noncircular steady states 
being elongated perpendicular to the field direction. 

In the following we will investigate the existence of 
steady states in the regime where z ~ 1. To this end, we 
expand the interior model in the small parameter 

8 = d+/R 

and look for first order perturbations of the steady state, 
i.e. 

R(9) = R + p(6)+O(S 2 ), 
c(r,9)=c + Cl +O(S 2 ). 
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Applying the quasistatic approximation, this leads to the 
following linear system for c\ and p: 

A Cl - r^Ci - 0, 



ci-d+C 1 c cos 9 = c° q -^(p + p"), 
= d t p = — («i — V ■ ri) 

2 r^ C l 



(14) 
(15) 



D(— x ■ ri - Vci • n ), (16) 

where no = (— cos 9, sin 0) is the normal of the circular 
steady state. 

With the ansatz (fTTj) for c\ , the steady state condition 
(fTBf is equivalent to 

O = d r /-^/cos0. 

Next, using (fT2|) and property (|A2|) of the Bessel func- 
tions /„ leads to the following simple recursion relation 
for the coefficients c n : 



Cn-lln-l + Cn+lln+l — C n {I n -l + ln+l) 



(17) 

where we have used the notation I m — / m (i?o/2£). For 
a solution that is symmetric under reflection at the x- 
axis (field direction) we have c„ = c_ n , which together 
with (|17p implies that c„ = cq. Therefore any symmetric 
solution of (fT4]) . (p~6j) is of the form 



ci(r, 9) = c exp(^-cos6») ^ I n (^)exp(in0) 



n— — oo 



= Co exp(| cos#), 

where in the second identity we have used (|A4j) . Now 
the boundary condition (|15p is used to fix the constant 
Co as follows: First note that (fTS"]) describes a driven har- 
monic oscillator in "time" 9, with the left hand side being 
the driving force. Next recall that a 2ir periodic solution 
p(9) of this oscillator exists provided the driving force is 
27r-periodic with vanishing n — 1 Fourier mode, where 
the latter condition means that the oscillator is not in 
resonance with the driving force. Using the Fourier ex- 
pansion of c\ (see (|A4[| ) this determines the constant Co 
to be 



Co 



Rp 



Thus, the steady state adatom concentration is given by 

a (r, 9) = d+ cq \ T *° exp( | cos 9) . 

Finally, the symmetric steady state shape (p(9) = p{—9)) 
is obtained as the solution of the ordinary differential 
equation |(T5j) as: 



p{9) = p {z) + Pn{ z ) cos n9 



iV ; n>2 



cosn6» . (18) 



The relative perturbation p/Ro is expressed as a function 
of the dimensionless parameters 8 = d+/Ro, z = i?o/C 
and T/Rq, which characterize the deviation from the TD 
limit, the strength of the electromigration force and the 
capillary effects, respectively. 




FIG. 2: First order perturbations of the circular steady state. 
Left: Steady state shape for different values of 8 = d+/ Ro <C 1 
with fixed z = R /£ = 5 and T/R = 0.05. Right: Different 
values of z with fixed 8z = 0.05. In both cases the perturbed 
shape Ro+ p — po is depicted, i.e. the dilation mode has been 
subtracted. 

The constant term po in (|18[) describes a uniform di- 
lation of the circle. Since T/Ro <C 1, po > 0, while 
Pn>2 < 0. In particular, the leading order deformation 
with n — 2 corresponds to an elongation perpendicular to 
the drift direction. This is illustrated in Fig. [2], where the 
dilation mode has been subtracted. Moreover, with in- 
creasing electromigration force (increasing z), the shapes 
start to become concave on the trailing side (recall that 
the islands are moving to the left). 

As has been pointed out above in Section HH the full 
(nonlinear) evolution is area conserving in the quasistatic 
limit. Since the dilation mode po is of the same order as 
the elongation mode p2, this property is generally vio- 
lated by the first order perturbation in 8. This suggests 
that the perturbative regime may be restricted to rather 
small elongations. To obtain the steady states of the 
full nonlinear model, numerical simulations of the time- 
dependent equations ©-(GO) have been performed. An 
adaptive finite element method is used, where the free 
boundary problem is discretized semi-implicitly using an 
operator splitting approach and two independent numer- 
ical grids for the adatom density c and for the island 
boundary, respectively. For the boundary evolution, a 
front tracking method is applied, for details see (33|. 

Starting with a circular shape, the void elongates until 
it reaches a steady state. A typical example of the time 
evolution is depicted in Fig. [31 Here the parameters are 
8 = 0.1, z = 5, which, as will be seen below, is already 
far away from the perturbative regime. We do not ob- 
serve any concave parts at the back side of the void as 
opposed to the perturbative steady state, see Fig.O This 
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FIG. 3: Simulation of the full nonlinear evolution of a circular 
vacancy island towards a steady state. For this example, S = 
0.1 and z — 5, which is well beyond the regime of validity of 
the first order perturbation theory. The void is moving from 
right to left and the snapshots are taken at t = 0, t = 1.5, 
i = 5, t = 50. The contours at later times are moved back to 
the right for better visibility. Space is measured in units of the 
radius Rq of the initial shape and time in units to = iio/|y |, 
where V is the drift velocity of a circular island as given in 
(|13|) . The simulation parameters are Ro = 100a, a 2 Cg q = 
10 -5 , F/Ro = 0.05. 



turns out to be true for all examples which we investi- 
gated. We also checked that the final steady state does 
not depend on the initial shape. For example, starting 
with a "bean-like" shape being the steady state shape as 
obtained from the perturbation theory also leads to the 
same convex steady state. Moreover, we have not seen 
any break up in the simulations even for cases where 
5 ~ 1. In all steady state shape similar to the 

one depicted in Fig. [3] is approached, where the deforma- 
tion increases with increasing 6 and where the curvature 
at the back side (right side) approaches zero. In Fig. 0] 
the steady states as obtained from the nonlinear evolu- 
tion are depicted for different values of 6. By comparing 




5 = 0.001 
5 = 0.010 
5 = 0.100 
5 = 1 .000 



FIG. 4: Steady state shape as approached by the nonlinear 
evolution for different values of S — d+/Ro <C 1 with fixed 
z = Ro/i = 5 and r/Ro = 0.05. 

with Fig. [5] (left) it is clear that the perturbative regime 



is limited to rather small values of 5, since already for 
8 = 0.01 the perturbative (non-convex) and the nonlin- 
ear (convex) steady state shapes differ considerably. 

In Fig. [5] we have investigated this quantitatively by 
comparing the deformation A 

A = (p(n/2) + p(-tt/2) - p(Q) - p(ir))/2Ro 

for the steady state obtained by the perturbation theory, 
i.e. given in ([Tg]) . and for the steady state as approached 
by simulations of the full nonlinear dynamics. The first 
order perturbation theory is seen to be quantitatively 
accurate only for 6 < 0.01. 




FIG. 5: Deformation A of the steady state shape versus rela- 
tive kinetic length <5 = d+/Ro, as obtained from perturbation 
theory and from the full model. 



IV. LINEAR STABILITY ANALYSIS 

Next we perform a linear stability analysis for the cir- 
cular steady state of the interior model in the TD limit. 
Thus we are looking for small perturbations of the form 

p(e,t) = p(6,0)exp(iut), 
ci(r, 9, t) = ci(r, 9, 0) exp(?w£). 

As in p^]) - p^|) we obtain - using the quasi-static approx- 
imation for c\ - the following linearized system for C\ and 
P, 

^CL-C 1 9xCl=0, (19) 

C i= c c q 4(p + P"), (20) 
dtp = a^D^-^x ■ no — Vci • no)- (21) 

In view of (|11|) we make the following ansatz 

p(6,0) =exp(f cos0)^p n exp(m0), z =~r- (22) 

71 
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Thus, the perturbation p{9, 0) is written as a series ex- 
pansion in the functions exp(| cos 9) exp(in0), instead of 
the usual Fourier modes. This choice has the advan- 
tage of simplifying the subsequent calculation. In par- 
ticular it will lead to a linear system, where each matrix 
row/column has only a small number of non zero entries. 
However, as we will discuss later, the non-orthogonality 
of the ansatz functions increases the non-normality of the 
matrix. Equation (|2ip now takes the form 



i(^2_j P n ex P(* n ^) = a2 D(d r f — t^/ cos#). 
Using (fT2|) , (|A2|) and matching coefficients we obtain 
a 2 D (, 



Ujjp n 



4i 



■\C n (ln-l + In + l) 



(23) 



where here and in the following we use the notation 
In = ^n(f)- We will finally use the linearized bound- 
ary condition (|20|) to express the coefficients c„ in (|23|) 
in terms of the p n 's. Inserting the ansatz (jlip for c\, and 
([22]) for p into ^ yields 



The right hand side of (|25|) represents the linearized time 
evolution of the island as a (infinite) matrix A acting 
on the coefficients p n - In real space it corresponds to 
an integro-differential operator, which is essentially non- 
local. Remarkably, the matrix depends on the system 
parameters only through the dimensionless electromigra- 
tion force z — i?o/£- I n particular, the capillary length 
r affects the time scale of the linear evolution, but not 
the stability of specific perturbations |l5|. For T — > 0, all 
eigenfrequencies uj n vanish, which implies that all pertur- 
bations become marginal. Indeed, it is easy to check that 
in the TD limit, for T = 0, any island shape translates 
rigidly at constant velocity under ([3])-([9]). 

The matrix A exhibits a symmetry, which originates 
from the invariance of the system under reflection at the 
x-axis (field direction). In terms of the coefficients p n 
this reflection is expressed as p n — > p~ n , and one readily 
verifies that this does leave (|25|) unchanged. Accordingly 
the eigenspace splits into two invariant subspaces with 
symmetric and antisymmetric eigenmodes characterized 
by 



Pn 
Pn 



P-n 
-P-n 



symmetric 
antisymmetric. 



f(Ro, 8) = c c ° q -^ J2{ 1 + Y {1 ~ COS20) ~ \ 
— inz sin 9 — n 2 ^jp n cxp(in9), 

which leads to (using (fT2]) for the left hand side) 

r 



C n In{ 2 ) — c eq 



— I (1 



z 



)Pn ~ -^(pn+2 + Pn-2) 



((2n-l)p„_ 1 -(27i + l)p n+1 )). (24) 



Inserting (|2"4"|) into (|2^|) leads after some tedious but 
straight forward calculations to the following eigenvalue 
problem for the coefficients {p n } and to: 

2 2 
Ap„ - (|(1 + y - n 2 )C n + y)p„ 

2 3 

+ {\n{ri + 2) + (2n + l)yC n - ^)pn+i 

2 3 

+ (^n(n - 2) + (-2n + l)yC„ - ^) p n ^ 



( 



with 
A 



^(2n + 3) + ^C„)/i„ , 

2 3 

+ (|g(2n - 3) -^C n ),,„_. 

+ + Pn-3), 



iuRl r _ In+l +Jn-l _ In+l 2rt 

a 2 c%DY , " ~ 27 n 7 n + 2 " 



(25) 



In both cases the eigenmodes are fully determined by only 
half of the coefficients (e.g. those with positive index), 
which allows us to reduce (|25|) to a semi-infinite system. 
By truncating it towards large n (cutoff towards small 
wavelengths) we arrive at a finite linear system which we 
solve numerically. 

Before we present the numerical results, we discuss 
translations and dilations, which are perturbations re- 
lated to the symmetry properties of the system. Symme- 
try under translations in the horizontal (x) and vertical 
(y) direction leads to two zero eigenmodes T x , T y - the 
infinitesimal horizontal and vertical translations - given 
by 

T x {9) = cos(0), T y (6) = sin(0). 

Indeed, inserting p — T x or p — T y into the linearized 
boundary condition ([2H|) leads to c% = and therefore 
(|2T1) yields dtp = 0. The horizontal translation T x be- 
longs to the symmetric eigenmodes, while the vertical 
translation T y belongs to the antisymmetric class. Next 
we consider a dilation T>, i.e., a constant initial pertur- 
bation 

V{9) = I. 

Inserting p{9) = T>(9) = 1 into the boundary condition 
([2"(If yields 



(-1 



c u r 

Rl 



which reflects the fact that one passes from one station- 
ary solution to another by increasing the radius of the 
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island and the concentration inside the island by a con- 
stant value. However, a dilation is not a zero eigenmode: 
Since we consider perturbations of a circular steady state, 
which has a steady state drift velocity depending on the 
radius, two circles with different radius are drifting apart. 
This leads to a linear increase of the perturbation. From 
(l2"Tj) . an initial perturbation p(9,0) = 1 has to grow ac- 
cording to 

a 2 c° YD , N a 2 c° YD , 
dtp = cos (0) = C 4—TJ6). 

In that sense a dilation T> generates a translation, and T> 
is a generalized eigenmode with eigenvalue zero according 
to A 2 T> ~ = 0. Thus (restricting to the symmetric 
case) the eigenvalue zero is two-fold degenerate, and has 
one (proper) eigenvector. Therefore the matrix A can not 
be diagonalized completely but contains a 2 x 2 Jordan 
block corresponding to the invariant subspace spanned by 
V and T x . Apart from that, the dilation doesn't play a 
role, because the time evolution preserves the area and we 
can therefore always restrict ourselves to perturbations 
which do not contain dilations. 




5 10 15 20 25 30 



z = 

FIG. 6: (Color online) Spectrum of the linearized theory as 
obtained by numerical solution of the eigenvalue problem 
Depicted are the eleven largest eigenvalues as a function of the 
dimensionless electromigration force z = Ro/£- Eigenvalues of 
symmetric modes are red/black, eigenvalues of antisymmetric 
modes blue/grey. Note that the eigenvalue zero is threefold 
degenerate. All negative eigenvalues come in pairs, consisting 
of a symmetric and an antisymmetric mode which become 
degenerate at z = 0. For the fluctuations appearing around 
z — 30 see the discussion at the end of Section IWl 

The numerically determined spectrum is presented in 
Fig. O Here the largest eigenvalues, for < z < 30 are 
depicted. The spectrum is purely real and apart from 
the predicted threefold degenerate zero eigenvalue it is 
strictly negative. For z — (no electromigration force), 
the right hand side of (f2l)|) becomes diagonal: 

Xp n = \n\(n 2 - l)p n , 



which allows to directly read off the eigenvalues A„ = 
\n\(n 2 — 1). Each eigenvalue A„ is twofold degenerate 
with eigenmodes given by the Fourier modes cos(n9), 
sm(n0). Since no field breaks the rotational symmetry, 
the symmetric and antisymmetric modes cos(n#), sin(n9) 
are connected by a rotation by 7r/2 and belong to the 
same eigenvalue A„ . In the presence of an electric current, 
i.e. z > 0, the rotational symmetry is broken and the de- 
generacy is removed, i.e., the eigenvalues split. Moreover 
increasing values of z lead to decreasing eigenvalues, i.e. 
larger islands or islands in the presence of a stronger field 
relax faster to the circular shape (as compared to the case 
without drift). For large values of z the eigenvalue prob- 
lem becomes numerically difficult to solve, leading to a 
noisy spectrum in Fig. [5] for values of z ~ 30. This will 
be discussed in more detail below. 

We now turn to the eigenmodes of the linearized time 
evolution. Figs. [7] and [H] show some examples of symmet- 
ric and antisymmetric eigenmodes with small index n for 
2 = 10. All of them reveal a typical undulated shape, 
where the number of nodes increases with the index and 
the amplitude is largest towards 9 = 0. Thus the modu- 
lation is more pronounced at the back side of the island 
(with respect to the drift motion) and, as shown in Fig.[9l 
this localization gets stronger with increasing z. This is 
a signature of the fact that the eigenvectors become more 
and more parallel (see below). The increasing localiza- 
tion of the eigenmodes at the back side of the island 
may be traced back to the convective, electromigration- 
induced flux in the adatom diffusion equation ([3]) which 
leads to the factor exp(| cos#) in l|22p. For large values 
of z, this factor strongly suppresses all contributions for 
values of 9, which are not close to 9 = 0. In particular, 
a perturbation that is front-back symmetric has to be a 
linear combination of many different eigenmodes. In that 
case, the eigenmodes with large index n will decay very 
fast leading to a shape being close to the steady state 
circular shape at the front side, while still having a large 
buckle at the back side. 

This behaviour will be investigated in more detail in 
the next section, when we consider the fully nonlinear 
evolution. For large values of z this will finally lead to a 
fingering instability at the back side of the island. 

Before turning to the nonlinear evolution, we comment 
on some observations connected with non-normal eigen- 
value problems. The matrix A in (|25p is highly non- 
normal for large values of z, which leads to a strong sen- 
sitivity of the spectrum with respect to small perturba- 
tions of the matrix entries. Small changes of the entries 
of the matrix may lead to structurally completely differ- 
ent spectra, a feature which can be consistently described 
within the theory of pseudospectra [3(|. The degree of 
non-normality depends on the choice of the basis. In our 
case the basis is not orthogonal (with respect to the usual 
L 2 scalar product) and therefore, although convenient for 
the analytical part, can give a distorted picture of the un- 
derlying geometry. We therefore checked the results by 
transforming the system (|2l)|) to the Fourier basis. Here 
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FIG. 7: Second, fourth and sixth symmetric eigenmode as 
function of the angle 9 (z = 10). 8 = 7r is the drift direction. 
The modes are normalized with respect to the L 2 -norm. 



FIG. 9: Symmetric eigenmode n = 3 for different values of z. 
The asymmetry with respect to front (6 = 7r) and back side 
(6* = 0) increases with increasing value of z. The modes are 
normalized with respect to the L 2 -norm. 




this and improve the condition numbers dramatically, of 
course by changing the spectrum beyond recognition. 

Finally we note that non-normality has also been dis- 
cussed [22|, [2H, [2t| [3(3] as a mechanism causing a transient 
growth of small perturbations even in the case of a lin- 
early stable system. This transient amplification can be 
large enough to drive the system out of the linear regime, 
this way leading to an instability. We thoroughly inves- 
tigated the transient behaviour of the linearized system 
and found that this scenario does not apply to our case. 



FIG. 8: Second, third and fourth eigenmode superimposed 
on a circle (grey background) for z — 10. Top (bottom) row 
shows symmetric (antisymmetric) modes. 



the non-normality is much less pronounced and we ob- 
tained the same spectrum as depicted in Fig. [5J In both 
cases, we have not been able to obtain the spectrum nu- 
merically for values z > 30 since, as can be seen from 
Fig. [6] the lines start to become noisy between z = 25 and 
30. The fluctuations are getting stronger very quickly, 
making it impossible to determine the spectrum reliably. 

A closer examination reveals that the condition num- 
bers for the eigenvalues grow exponentially with z and 
therefore the accuracy of the results gets lost very quickly 
at some point. This is again a consequence of the non- 
normality of the problem. The condition number of a 
given eigenvalue is large when the corresponding left and 
right eigenvectors are almost orthogonal [34]]. This is in- 
timately connected with the non-orthogonality between 
different (right) eigenvectors. In our case the eigenvec- 
tors become almost parallel for large values of z. Even 
a slight perturbation of the matrix entries will destroy 



V. NONLINEAR EVOLUTION 

To further investigate the nonlinear evolution, we 
have performed numerical simulations of (O)-© using 
an adaptive finite element method, as described in Sec- 
tion [TTTJ We note that adaptive mesh refinement during 
the time evolution in regions with high curvature turned 
out to be crucial for accurate simulations in the case when 
small fingers appear in the shape. 

We probe the nonlinear dynamics with a very strong 
electromigration force, i.e. choosing z = Ro/£ — 25. In 
Fig. [10] a simulation of a vacancy island (interior model) 
in the TD regime is depicted. Here an ellipse with aspect 
ratio 1.2 has been taken as the initial shape. We observe 
a fingering instability at the back side of the vacancy 
island, which finally leads to a pinch-off. Note that this 
behaviour has not been predicted by the linear stability 
theory. 

To gain a better understanding of the nonlinear dy- 
namics, the initial evolution of an ellipsoidal island has 
been investigated in more detail. For this special geome- 
try, one may use elliptic coordinates and an expansion of 
the adatom concentration in terms of Mathieu functions 
to arrive at an analytic expression for the normal veloc- 
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FIG. 10: Nonlinear dynamics in the interior model with 
strong electromigration force (z = 25). The initial shape is an 
ellipse with aspect ratio 1.2. A fingering instability appears 
at the trailing end, which leads to a pinch-off. The contours 
are shifted upwards for better visibility and also shifted to 
the right with Ax = t\V\, where V is the drift velocity of a 
circular island, see (|13|l . Space is measured in units of the 
lattice spacing a and time in units to = 7?o/|V|. The simula- 
tion parameters are Ro — 100a, a 2 c° q = 10 -5 , T/Ro = 0.02. 
The breaking of the up-down symmetry in the right column 
is due to numerical noise. 



ity, see Appendix [Bl (|B2[) . However, as it turned out, 
a reliable evaluation of the series expansion is possible 
only for moderate values of z < 5 and aspect ratios < 2. 
For parameters in this regime, the normal velocity of an 
ellipsoidal vacancy island - elongated in the horizontal 
direction - in the center of mass system, i.e. after sub- 
tracting the (fairly large) drift velocity, is compared with 
the corresponding results of the finite element simulation 
in Fig. HU 

The two solutions are in very good agreement. More- 
over, one realizes that the dynamics is fastest at the front 
part of the island and tends to relax the front half of the 
ellipse towards a circle, since the boundary moves in- 
wards at = n ( front ) but outwards in nearby regions 
with 9 ps 7r ± 7r/4. Since the evolution at the back half 
of the island is slower and always inwards, we expect the 
dynamics to lead to an egg like shape, which is verified in 
the simulations, see Fig.[T2] Moreover, the asymmetry of 
the dynamics becomes more pronounced with increasing 
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FIG. 11: Normal velocity of an ellipsoidal vacancy island 
(elongated in the horizontal direction) as a function of the 
normal angle in the center of mass system, i.e., after sub- 
tracting the drift velocity of the island. The normal points 
inwards, i.e. a positive velocity corresponds to local shrinking, 
and the angle 9 — n corresponds to the drift direction. The 
normal velocity taken from the numerical simulations by aver- 
aging over the first 1000 time steps is in very good agreement 
with the results obtained by evaluating the analytic solution 
using elliptic coordinates. 



values of z. We expect this to be a possible reason for 
the onset of a fingering instability, once a critical value 
of the curvature at the back end is reached. 

To roughly locate the threshold for the onset of the 
instability we systematically varied the electromigration 
strength z and the amplitude of the initial perturbation 
given as the aspect ratio of the initially horizontally elon- 
gated circle. The results of the simulations are presented 
in Fig. [T2J As can be seen, the instability sets in around 
z ~ 10. 

Although we are so far lacking a satisfactory analytic 
understanding of the underlying mechanism, we would 
like to present some intuitive reasoning for the onset of 
an instability at large values of z. First recall that the 
mass flux at the boundary of the vacancy island in the 
normal direction n (pointing inwards) is the sum of the 
diffusive flux 
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and the flux caused by electromigration 
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where, as usual, x denotes the unit vector in ir-direction. 
Next we consider the two simplified cases of (a) j e = 0, 
i.e. no electric field and (b) jd = 0, i.e, no diffusion, to 
argue that jd is stabilizing while j e is destabilizing: 

(a) No electric field: An outward bump of the bound- 
ary leads to a local maximum of the curvature and 
therefore a local minimum of the adatom density 
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FIG. 12: Numerical simulation of the evolution of an ellip- 
soidal vacancy island with aspect ratio 2.0 and area tvRq. 
The upper figure depicts the initial stage of the evolution: 
the left/right asymmetry increases with increasing value of 
z = Ro/£, leading to the formation of an egg-like shape. Bot- 
tom: For large values of z, the asymmetric evolution leads to 
a fingering instability at the back side of the island. Space is 
measured in units of the lattice spacing a and time in units 
to — Rq/\V\, where V is the drift velocity of a circular is- 
land as given in (|13p for the case z — 10. The simulation 
parameters are Ro = 100a, a 2 c[! q = 10~ 5 , F/Rq = 0.05. 
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FIG. 13: Threshold for the onset of the instability: Simu- 
lations for different values of z = Ro/£ and aspect ratios of 
the initial ellipsoid, where the elongation is in field direction, 
have been performed. The instability sets in at z — 9. For 
z > 11 also the smallest aspect ratio leads to an instability, 
i.e. a fingering at the back side of the island. 



c eq along the boundary. Due to the maximum prin- 
ciple, this will also be a local minimum of c inside 
the island, leading to a diffusive flux towards the 
boundary, i.e., the bump will be filled. Thus the 
diffusive flux is stabilizing. 

(b) No diffusion: Consider an outward bump of the 
boundary at the trailing (back) island edge. As 
in case (a), this leads to an increased curvature 
and therefore a decreased adatom density c cq as 
compared to the front side. The flux j e caused by 
electromigration is proportional to c cq and there- 
fore decreased at the back side, i.e. this part of 
the boundary is now drifting more slowly than the 
center of mass and the bump is growing. Thus the 
electromigration flux destabilizes the trailing edge 
and stabilizes the leading edge of the island. 

Finally consider the general case of an outward bump 
of the trailing boundary in the presence of diffusion 
and electromigration. Fixing all parameters, except the 
strength of the electric field, the diffusive flux de- 
pends essentially on the geometry of the bump, whereas 
the decrease of j e is proportional to £ . Thus we may 
expect that the bump will grow if the electric field is large 
enough. 

The argument also elucidates the role of the capillarity 
parameter T in this problem, which is mainly to trans- 
late variations of the boundary curvature into variations 
of the adatom concentration in the interior of the island. 
The latter in turn underlies both the stabilizing effect of 
the diffusive current and the destabilizing effect of elec- 
tromigration. In contrast to the instabilities in the PD 
regime, which can be described in terms of a competition 
between capillarity and electromigration [22|, HE HtJ , here 
capillarity plays a largely neutral part, which explains 
(in hindsight) why the linear stability properties are in- 
dependent of T (see Section IIV[) . 



VI. CONCLUSIONS 

In this paper, we have presented a detailed study of a 
continuum model for the electromigration-driven shape 
evolution of single-layer vacancy islands mediated by in- 
ternal terrace diffusion. Significant shape deformations, 
such as the elongation transverse to the field direction de- 
scribed in Sect ion [HI] or the pinch-off instability discussed 
in Sect El were found to require dimensionless electromi- 
gration forces z = Rq/S, significantly larger than unity. 
Unfortunately this implies that these phenomena will be 
difficult to realize experimentally, at least for the surfaces 
commonly used in this context. To give an example, the 
maximal electromigration bias that can be achieved on 
the Cu(100) surface has been estimated [35[ to be on the 
order of -Ebias ~ 10~ 5 eV for a diffusion hop between 
nearest neighbor sites, corresponding to a characteristic 
length scale £/a = k^T/E^^ 2.5 x 10 3 . An island 
with z = 10 would thus contain about 2 x 10 9 atoms, 
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which is four orders of magnitude larger than the size at 
which strong shape deformations and electromigration- 
induced oscillatory dynamics have been predicted in the 
PD regime [H. 

From the broader perspective of the theory of moving 
boundary value problems, our results are of interest be- 
cause they add another example to the list of cases in 
which the standard tool of linear stability analysis fails 
to correctly predict the stability properties of the full 
nonlinear dynamics. The behavior of the interior model 
in the TD limit is similar in many respects to the void 
mi grat ion problem in the PD regime which was studied 
in [22L [H| ■ In both cases the circular solution is lin- 
early stable for arbitrary values of the electromigration 
force, but a nonlinear instability occurs when the sys- 
tem is subjected to finite amplitude perturbations, and 
the threshold for the nonlinear instability decreases with 
increasing force. Moreover, as in the case of void migra- 
tion, a distorted island either relaxes back to the circular 
shape or evolves towards pinch-off. This suggests (as has 
been proven for voids [251 ]) that no non-circular station- 
ary states exist. However, in contrast to the void migra- 
tion problem, in the present study we found no evidence 
for transient amplification of linear perturbations related 
to the non-normality of the eigenvalue problem. 

In this context it is worth mentioning the problem of 
two-dimensional ionization fronts, which shares some of 
the features of both void and island migration [37|]. In 
this system a closed curve representing the ionized region 
(a "streamer" ) evolves in response to a Laplacian poten- 
tial in the exterior domain. As in the void migration 
problem, a constant potential gradient far away from the 
streamer represents the driving electric field. However, 
the boundary condition at the front is similar to that 
employed in the present work, (|5I6I9[) with c cq = const., 
i.e. r = 0, and the kinetic length corresponds to the 
width of the ionization front. Again, circles translating 
at constant velocity are solutions of the problem. In the 
special case where the kinetic length equals the radius 
of the streamer the linearized dynamics can be solved 
exactly, and it is found that the circle is always stable. 
Standard linear stability analysis nevertheless fails, be- 
cause smooth initial perturbations do not decay expo- 
nentially [38| . Further exploration of the relationship be- 
tween these three boundary value problems seems like a 
promising direction for future research. 

APPENDIX A: BESSEL FUNCTIONS 

We summarize some properties of the modified Bessel 
functions of imaginary argument /„. For integer n the 
functions /„ are symmetric with respect to the index 

J_„(r) = I n (r). (Al) 

The derivative is 

| : / n (r) = i(7„_ 1 (r)+J n+ i(r)). (A2) 



There is the recursion relation 

r/ n _i(r) - rl n+1 (r) = 2nl„(r), (A3) 

and the generating function is 

00 

exp(rcos6*)= I n (r) exp(in6). (A4) 



APPENDIX B: NORMAL VELOCITY OF AN 
ELLIPSOIDAL ISLAND 

To calculate the concentration profile inside an ellip- 
soidal island, we introduce elliptic coordinates 

x = acosh(u) cos(u>), 
y = asinh(u) sin(w), 

where u and w are the radial and angular coordinates, 
respectively. The line u = uq = const, is an ellipse with 
aspect ratio tanh(uo) and the parameter a determines 
the size of this ellipse. Curvature and normal vector are 

k = 5- cosh(uo) sinh(uo), 

_ _ 1 / sinh(uo) cos(w) \ 
U ~ ^fg V cosh(w )sin(w) J ' 

with g — cosh 2 (uo) — cos 2 (w). The Helmholtz equation 
now reads 

(dl + d 2 Jf = a 2 (sinh 2 ( W ) + sin 2 ( W ))/ 

and separation f(u,w) — f u (u)f w (w) leads to 

fw ~ » 2 sin 2 (w)f w = Xf w , 
~f: + a 2 sinh 2 ( U )/ M = Xf u . 

The solutions of the two equations are related by 
f w (ix) = f u (x) . Requiring periodicity of f w determines a 
discrete set of values for the separation parameter A. The 
corresponding solutions are the Mathieu functions of the 
first kind ce n and se n [39[ . Since the ellipse is symmetric 
with respect to the field direction we only need the ce n 
which are even functions of the angular variable w. The 
general solution is then 

c(w, w) — exp ( ^- cosh(w) cos(u;) J b n ce n (w)ce n (iu). 

(Bl) 

The coefficients b n are determined by the boundary con- 
dition (lU). The Mathieu functions are orthogonal and 
normalized according to L"ceJ(i) dx = n. Thus, the 
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b n can be found in a way similar to Fourier coefficients 
via the integrals 



K. = 



2jv 



7rce„(mo) Jo 
x (1 — Tk) dw. 



ce„(w)exp ( — — cosh(u ) cos(ui) ) x 



These integrals are solved numerically. Finally, using the 
general solution (IB1[) to calculate the flux to the bound- 
ary, the normal velocity v of the island edge is obtained 
as 



exp ( — cosh(uo) cos(u/ 



E 



b n ce n (w)x 



_ a 2 D 

x ( — Im(ce^(mo)) + — sinh(ito) cos(ui)ce„(iuo 



a 



(B2) 
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